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CONVECTIVE INSTABILITY OF A POLYTROPIC LAYER IN A 

MAGNETIC FIELD AND A CORIOLIS FORCE 

By Kin-Wai Yeung* 

Manned Spacecraft Center 

SUMMARY 


The onset of convection induced by gravity is considered in a horizontal, infinite, 
polytropic fluid layer, heated from below, under a uniform vertical magnetic field and 
a coriolis acceleration resulting from a uniform rotation. The viscosity and thermal 
and electrical conductivities are included in the analysis. The principle of the ex- 
change of stabilities is assumed valid. The numerical method is used to compute the 
critical Rayleigh number that defines the onset of the instability. It is found that the 
stability of the system depends on the layer thickness. 

When the rotation rate is small, stability increases as the magnetic field in- 
creases. However, when the rotation rate is relatively large, complex interactions 
occur between the rotation and the magnetic field. The critical Rayleigh number 
decreases with increasing magnetic field until a minimum value is attained. From 
that point, the critical Rayleigh number increases as the magnetic field increases. 

For a small magnetic field, rotation is effective in stabilizing the system. How- 
ever, if the magnetic field is large, the stabilizing effect achieved by increasing the 
rotation rate becomes less prominent. 


INTRODUCTION 


The analysis of the thermal instability induced by gravity in an infinite layer of 
fluid medium is the subject of numerous examinations, because of the potential appli- 
cations to both geophysical and astrophysical phenomena. Analyses using models 
without rotation have been studied extensively (for example, refs. 1 to 4). However, 
in the application to a stellar atmosphere, a coriolis force should be included in the 
study. Therefore, many investigations that include the rotation have been made. An 
incompressible medium, subjected to a coriolis force and a magnetic field simultane- 
ously, has been examined by Chandrasekhar (ref. 5) and Talwar (ref. 6), who also 
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studied the compressible model but considered the medium as an ideal fluid having 
infinite electrical conductivity (ref. 7). A similar model was analyzed in other studies 
(refs. 8 and 9) that included compressibility but ignored the electrical nature of the 
medium. 

To provide a more realistic picture of thermal instability from the astrophysical 
standpoint, the finite transport coefficients, as well as the compressibility of the me- 
dium, have been included in the investigation discussed in this paper. 


SYMBOLS 


A.,B. 

l l 


a 

b 

C 


c 

P 

c 

V 

D 

E 

e 


g 

H 


H 


H 


H 


1 ~'F' 

J 


coefficients used in equation (45) 

horizontal wave number defined by equation (37) 
layer thickness 

coefficient defined by equation (15) 
specific heat capacity, constant pressure 
specific heat capacity, constant volume 

d 

operator defined as ^ 
electric field 
exponential 
gravitational field 
magnetic field 

x-component of the perturbation magnetic field 
y-component of the perturbation magnetic field 
z-component of the perturbation magnetic field 

current density 
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K thermal conductivity 

k horizontal wave number, defined by equation (33) 

k wave number in the x- direction 

x 

k wave number in the y-direction 

y 

m parameter defined by equation (14) 

P integration constant 

p pressure 

Q Hartmann number 

Qj parameter defined by equation (63) 

R Rayleigh number 

R c critical Rayleigh number 

R* universal gas constant 

T dimensionless thickness 

Taylor number defined by equation (62) 

Tj parameter defined by equation (64) 

r parameter defined by equation (40) 

t time 

u,v,W x-, y-, and z-components of the perturbation velocity 
W component defined in equation (27) 

x,y,z rotating coordinates 

0 thermal gradient 

T polytropic index 

e variable defined in equation (42) 
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? dimensionless variable defined by equation (35) 
77 electrical resistance, defined by equation (23) 

0 dimensionless temperature 

jj. material viscosity 

p. magnetic permeability 

p density 

a electrical conductivity 

t stress tensor 

£2 angular velocity 

co amplitude of lo , defined in equation (29) 

z 

co z- component of y Xu 

z 

( ) vector 

( )' perturbation quantity 

( ) tensor 

A 

( ) unit vector 

( ) amplitude of the perturbation quantity 

Subscripts: 

0 equilibrium condition 

1 index, i = 1, 2, 3, ... 

Operator: 
v del- operator 
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GOVERNING EQUATIONS 
OF THE PROBLEM 


x 


The system considered in this paper 
is a free infinite layer of thickness b, ro- 
tating with a constant angular velocity £2 
about the vertical axis. The configuration 
is shown in figure 1. Both the applied 

magnetic field Hq and the gravitational 

field g are parallel to the z-axis and 
are considered to be constant. 

In rotating coordinates x, y, 
and z, the governing hydrodynamic equa- 
tions are 


cu ■e T t 2_i 

P at p ( u ' v ' u " ^ e J X H = " vP + V-V u 

+ PS + -g |iv(v • u) + 2pu X n (1) 



Figure 1.- Configuration of the system 
with a free infinite layer of thickness 

b, rotating with angular velocity £2 
about the z-axis. 


|f + v • (pu) = 0 


( 2 ) 


pc y + (u • v)tJ= Ky 2 T - Pv • u + -^f- + (t • v) • U + p.X7 • (u • vu) 


I d(v 




( 3 ) 


and the equation of state is 


p = R*pT 


( 4 ) 


The equations governing the electromagnetic conditions are 


( 6 ) 


V 


X E = -fi 


9H 
e at 


and 


v • H = 6 



( 7 ) 


t 

(8) 


In equations (1) to (8), p is the material viscosity, p e is the magnetic permeability, 

and t is the stress tensor, which is a function of derivatives of u. The electromag- 
netic units are used in this paper. 

From the static equilibrium 


P 


0 " 


Pz 


m+1 


(9) 


„ Pm 

(“>) 

and 

T o = V dD 

where P is an integration constant, /3 Q is the applied thermal gradient, and the 
origin of the coordinates is chosen at the level where Tq = 0. Variables with a sub- 
script zero represent the equilibrium conditions. Equations (9) to (11) are indicative 
that 


P 


0 



(12) 
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■ ill iiiii ■mfeiifliiiMllimiMM niiillll i i 



with 


T = — ^ - = polytropic index 


m “ /3 q R* ' 1 


C = Pi 


l+(l/m) 


Hence, the equilibrium state follows the polytropic relation. 


THE PERTURBATION STATE 


Consider a slightly perturbed state in which the variables can be expressed as 


P = P 0 + P' 


P = P 0 + P 


T = T Q + T' 


u = u' = (u, v, W) 


E = E' 


J = J’ 


H = H 0 + H’=( H i, H ;. H o + H;) 
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The primed terms are the perturbation quantities that have negligible products. 
Assuming that the marginal stability is valid (that is, the a/at terms vanish), the 
governing equations can be combined to form the following equations. 


and 






4tT \ H 0 3z/ H ' = M p ' ¥zj + ^“u’ + P’g 


Mv(v • u’) + 2p 0 O(ve x - ue J 


v • u' = 


— W 

Z 


(17) 

(18) 


% c W = K v 2 T’ 


(19) 


V'*’-( H o^ 4 ’ = w 2 ff 


( 20 ) 


P’ = 


i_ L _ pz 


m 


R*T, 


P’ - 


( 21 ) 


/ 3 = ^n- 


g 


0 c 


( 22 ) 


47 ll± _o 


(23) 


When standard eliminations by means of curl and divergence operations are 
used on equation (17), the results can be combined with equations (18) to (21) to give 


M e H o' 4){a t £ w _ „ 


47777 


m + 1 dz \z dz 


2 12 z a 

v + 3 V 1 ~ nT+T v dz. 


+ V 


2 a 


3z m + 1 


V 


Pz m 

w -^-v 

p 0 


® l 1 - iTTl 5?) V, = 0 <*«> 
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(25) 


V 



MeH ° CO + 2- JL 

4nnrj dz 2 z n dz 



= 0 


and 

/3p c W = K V 2 T’ (26) 

^0 p 


2 2 2 2 2 

where a> is the z-component of v Xu 1 , and Vi - 9 /3x + 3 /3y . 

Z 1 

Because the coefficients in equations (24) to (26) are not functions of x and y, 
the solutions are assumed to be in the form 


W=W(z)e‘( k x x + V) 

(27) 

T'» 6 (z)e I ( k x x+ V) 

(28) 

to = ®(z)e'( k x x + k y y ) 

(29) 


where k x and are the wave numbers in the x- and y-directions, respectively. 

By substituting equations (27) to (29) into equations (24) to (26), the following equations 
are obtained. 



PR* 
M g 


/ i\2 , 2 m-1 

(m + 1) k z g 


20 

M 


(^ 



p 0 " = ° 


(30) 
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d 2 2\ _ M e H 0 d 2 _ „S2 d 

2 - kj “'-BiS- “ +2 ^ d^(^ 


Vdz 


^(4 - ^ M - 0 (3i > 


and 


/3c„ ^ W = k/- 4 - k 2 ) 0 


'P R*/3 f 


Vdz 


(32) 


with 


2 2 2 

k = k + k = horizontal wave number 

x y 


(33) 


By use of equations (30) to (32), the following equation governing e can be 
formed. 


l)m 2 

— — a 


| D(ljD 2 [|D(0](j- (2a 2 + Q)Dji[i D (0][ + ^ 

1 FL2 2'| 2 2'\ /m\ 1 (m + 

||D -a) + (D " a /D^j +3 “ ^ 

. (2L±i) ( D 2 - a 2 ) (d + ”) - q(d - + d) 

+ r ’D | d(d 2 - a 2 )(c m j | [? _m (D 2 - a 2 )]© 

= -?0" (2m " 1)Ra2 | D [|D 2 (iD)] - (2a 2 + q)d^d) + - 


(34) 


This equation has been normalized by the layer thickness b, and the dimensionless 
quantities are defined as 



( 35 ) 


10 



(36) 


and 


D = 


‘(5) " dC 


2 2 2 

a = k b = horizontal wave number 


(37) 


2^., n / ., >3 ,2 m +3 

„ . 2m-l P + b 

R ^0 „ 2 

K-ig 


= Rayleigh number 


(38) 


TT 2, 2 
|a e H 0 b 

Q = = Hartmann number 

4*117) 


(39) 


r = b 4 1— ]2 f Pbm 


u/ v/3 0 R *y 


(40) 


THIN-LAYER ANALYSIS 


Because equation (34) is a lOth-order ordinary differential equation with variable 
coefficients, an analytical solution is difficult to obtain. For simplification, the anal- 
ysis will be restricted to the thin-layer case. 

Consider a layer that has a small thickness b and is located at a large distance 
Zq from the origin such that 


T = 




(41) 


The parameter T is defined as the dimensionless thickness. The variable £ can 
now be written as 


S = ? 0 +e 


(42) 


11 



where e varies from zero to one. Consequently, quantities in the form of l/£ 
be expressed as 


S 


n "TT fl - n -J- 

n „ n | 


and 


d_ _ d 
d? " de 


With this approximation, equation (34) can be written as 


,10 .9 ,8 ,7 ,6 

A ® + A 1-6 + a 1-i + A 1^ + A 1J + A de 

-* r\ ■< r\ u” n ' ^ O O '* Ft - “r -Tl/j £> r -tt.,- 


10 , 10 9 , 9 8,8 7 7 7 Tn 6 J T ~5 , 5 

de de de de de de 


+ ( a 4 + Ra b 4 ) ^4 + ( a 3 + Ra b 3 ) ^3 


2 „ \ de 


+ (A 2 + Ra2B 2 )fj + ( A i + Ra “ B i;af 


+ I Aq + Ra B 0 Je = o 


where 


A„ = -a 10 + 5a 10 J- + ^-p (l2a 8 + 4Qa 6 ) 

0 ? 0 ?o 2 

+ e ( 34 a 8 + 28Qa 6 ) + -V (236a 6 + 288Qa 4 + 84Q 2 a 2 ) 
Cn ?n 


-^(208a 6 + 248Qa 4 + 72Q 2 a 2 ) - (52a 4 + 36a 2 Q) 


can 


(43) 


(44) 


(45) 


(46) 
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A, = - -J- (25a 8 + 13Qa 6 + a 4 r) + (74a 8 + 40Qa 6 ) 

1 ? 0 S„ 2 


- — L (l!2a 6 + 130Qa 4 + 38Q 2 a 2 ) 


+ -^j (660a 6 + 624Qa 4 + 180Q 2 a 2 ) 


(858a 4 + 510Qa 2 ) - (444a 4 + 252Qa 2 ) - 


21 72a 


A 2 = (5a 8 + 2Qa 6 + a 4 f)- (25a 8 + 10Qa 6 + a 4 t) 


-^2 (l70a 6 + 118Qa 4 + 30Q 2 a 2 ) + ~ (24a 6 + 10Qa 4 + 12Q 2 a 2 ) 


(l354a 4 + 958Qa 2 + 84Q 2 ) + (l794a 


4 + 1166Qa 2 + 72Q 2 ) 


(712a 2 + 36 q)- 


A„ (lOOa 6 + 64Qa 4 + 13Q 2 a 2 + 2a 2 t ) 


(298a 6 + 194Qa 4 + 40Q 2 a 2 ) 


+ (l50a 4 + 146Qa 2 + 38Q 2 ) - (l 


4 (,1516a 4 + 1124Qa 2 + 180Q 2 ) 


\ (2112a 2 + 510 q) + (l488a 2 + 252 q) + 
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A 4 = -(lOa 6 + 6a 4 Q + Q 2 a 2 + 2a 2 <C ) 

+ J- (50a 6 + 30Qa 4 + 5Q 2 a 2 + 2a 2 T- ) 
^0 

+ ^2 ( 438a4 + 262Qa 2 + 30Q 2 ) 


5, 


(294a 4 + 178Qa 2 + 12Q 2 ) 


+ (l672a 2 + 670 q) - ~ (2672a 2 + 918 q) - +216 — 

to ^ to 


A 5 = -f (l50a 4 + 89Qa 2 + 13Q 2 + r ) 


+ (450a 4 + 270Qa 2 + 40Q 2 ) + — (30a 2 - 16 q) 

+ (l 048a 2 + 500Q) + i 2 ^. _ 1044 e 

Cn C„ 5 ° 


A = (lOa 4 + 6Qa 2 + Q 2 + v)- y- (50a 4 + 30Qa 2 + 5Q 2 + r 

5 0 


A ? = (lOOa 2 + 38 q)~ - (302a 2 + 116Q 


e 70 192e 


■8 = -( 5 ‘ 2 


A n = -l5a“ + 2 q)+ 5(5a 2 + 2 


oM + 


134 144e 

? " 2 - 3 
u ^0 ^0 


Ag =--1^ + 76 

? 0 ^ 


(50) 


(51) 

(52) 

(53) 

(54) 

(55) 
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( 56 ) 


and 


A - i _ 

10 _1 TZ 


4 4 p 

B Q = a - 2a 


0 


(57) 


Bj_ = 2 ( 2 a 2 + Q )* 




(58) 


B„ = -(2a 2 + q) + |l-(2a 2 + q) + — 
^0 


(59) 


B - ± + _?L 

^3 " " + „ 2 


'0 ? 


0 


(60) 


R - 1 2e 
b 4 " 1 “ r 
^0 


(61) 


Note that 7“ in equation (34) is now replaced by V which is defined as 


4/2m 2 | Pz 0 


m\ 


u. 


r = b ( — ) I j - Taylor number 


(62) 


For convenience in comparing the results with those obtained by Chandrasekhar 
(ref. 5) for the incompressible model 


Q = — 

^1 2 

TT 


(63) 
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and 


r 



(64) 


will be used in the computations, for which m is set equal to 2. 


BOUNDARY CONDITIONS 


Equation (45) must be solved with appropriate boundary conditions. Assuming 
that the perturbations vanish at the boundaries, the following equation is derived. 


W = e = 0 at e = 0, 1 


(65) 


Because only the free layer is considered, no boundary stresses exist at the surfaces. 
This condition leads to 


d 2 w + ™ W = 0 

(66) 

DH = 0 

(67) 

DoS = 0 

(68) 

D 55 = 0 at e = 0, 1 

(69) 


The adjoining vacuum acts as a perfect insulator, thus requiring that 


E = 0 at e = 0, 1 
z 


(70) 


Because of the variable coefficients of the equations, the other boundary condi- 
tions for e are difficult to obtain. To simplify the boundary conditions, all terms 

that contain lA Q n , where n >1, will be neglected. Only the thin- layer analysis is 
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considered; therefore, such an approximation is justified, because l/Cg 0 is much 

smaller than unity. The appendix contains a discussion of the use of approximate 
boundary conditions for the thin- layer analysis. 

The following two points must be clarified: (1) In deriving the governing equa- 

tion (eq. (45)), all l/? n terms are approximated with their linear forms obtained 
from the binomial expansion, as shown by equation (43). By using this method, the 

compressibility of the medium is retained. If all l/£ n terms are allowed to vanish, 
the problem is reduced simply to that of the incompressible model used by Chandrasek- 
har (ref. 5). (2) For the boundary conditions, the l/? n terms for the incompressible 

model are used. The approximation greatly simplifies the problem and is justifiable 
for the thin-layer case. Because >> 1, as defined by equation (41), £ = con- 

stant, and Pq will be treated as constant in deriving the boundary conditions. Then, 
equation (30) can be normalized by b and can be written in the form 


(d 2 - a 2 ) 2 - QD 2 


7 F 7 PR*/ ,\2 2, -m+3 „ m+1 
W + ——{m + 1) a b C 0 t 


20 Pb" m+3 

n R *^ 0 


? 0 m Dw = 0 


(71) 


Equation (66) is then reduced to 

D 2 w = 0 at e = 0, 1 (72) 

By use of equation (65), equation (32) can be shown as 

D 2 e =0 at e = 0, 1 (73) 

By differentiating equation (32) twice and combining the result with equations (72) and 
(73) 

D 4 q = 0 at e = 0, 1 (74) 
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When equations (65), (68), and (72) are substituted into equation (71), the result is 



0 at e = 0, 1 


(75) 


When the fourth differentiation of equation (32) is combined with equations (74) 
and (75) 


D 6 e 


= 0 at 6=0, 1 


(76) 


If equation (71) is differentiated twice and equations (69), (72), and (75) are used, the 
result is 


d 6 w = 0 


at e = 0, 1 


(77) 


The substitution of equations (76) and (77) into the sixth differentiation of equation (32) 
gives 


D 8 e = 0 at e = 0, 1 


(78) 


Hence, the 10 approximate boundary conditions for q are 


9=0 

> 




D 2 0 = 0 





D 4 0 = 0 



at e = 0, 1 

(79) 

o 

11 

CD 

Q 






and 

D 8 0 = 0 


These boundary conditions, together with equation (45), form an eigenvalue prob- 
lem, having and as the parameters and the Rayleigh number R as the charac- 
teristic value. Because R indicates the stability of the system, the minimum value 
R £ , with respect to the horizontal wave number, defines the onset of the thermal 

convection. 
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RESULTS 


Equation (45) can be reduced to the 
equation obtained by Chandrasekhar 
(ref. 5) by setting T = 1 /Cq = 0. Hence, 

the results obtained for the incompressible 
case correspond to the thin- layer case 
with infinitesimal thickness. 

As shown in figure 2, the critical 
Rayleigh number R c decreases as the 

dimensionless thickness T increases for 
a given value of Qj and Therefore, 

in a compressible model, stability de- 
pends on the layer thickness. Because it 
is assumed that the perturbations vanish 
at the boundary, a thicker layer therefore 
allows more freedom for the disturbance 
than the thinner layer. 

As shown in table I. for a given 
thickness and a given small magnetic 
field Qj. the rotation rate increases 


6.4 X 10^ 



Figure 2. - Critical Rayleigh number as a 
function of dimensionless thickness T. 


TABLE I. - CRITICAL RAYLEIGH NUMBER AS A FUNCTION OF Cj FOR Qj = 10, 


Qj = 500, Qj * 10 000, AND T = 0. 06 


®1 


R c 

a 

10 

1 

2 351 

3.60 


50 

3 611 

4. 00 


100 

4 848 

4. 30 


500 

12 802 

6.76 


1000 

20 009 

8.10 

500 

1 

55 350 

7.50 


50 

55 421 

7. 50 


100 

55 493 

7.50 


500 

56 337 

7.50 


1000 

57 054 

7. 50 

10 000 

1 

944 576 

12.6 


50 

949 756 

12.6 


100 

949 765 

12.6 


500 

949 837 

12.6 


1000 

949 926 

12.6 
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the stability, because R increases rapidly as increases. The size of the cell 

also increases accordingly, as indicated by an increasing a. However, for a large 
value of Qj, a large increase of only increases the stability to a slight extent, and 

the size of the cell remains the same. An increased rotation rate can promote stabil- 
ity because the vortexes generated by the rotation increase the tendency of the fluid 
particles to remain at the same altitude and to resist being moved to a higher level by 
convection. 

When the magnetic field is very large, the fluid particles are strongly bound to 
gyrate about the magnetic lines and have less drift tendency. The circulating motion 
of the fluid element is caused primarily by the Larmor gyration rather than by the 
applied rotation. Consequently, the effect of the applied rotation upon the stability be- 
comes less dominant. 

The combination of the magnetic field and the rotation rate helps to delay the on- 
set of the thermal convection (table I). However, this statement is not exactly true 
when ^ i is very large, as indicated in table n. For^j » 100 000 and for 

TABLE H. - CRITICAL RAYLEIGH NUMBER AS A FUNCTION OF Qj 
FOR = 100 000 AND T = 0. 06 


*1 

a 

R c 

20 

19.2 

399 111 

30 

19.1 

398 336 

40 

19.0, 3.6 

397 460, 876 700 

50 

18.8, 3.5 

396 589, 719 920 

60 

)— 1 

00 

00 

CO 

395 710, 611 410 

80 

18.8, 3.3 

393 982, 472 281 

100 

18.6, 3.3 

392 156, 387 757 

200 

3.4 

222 611 

500 

4.2 

154 566 

1 000 

5. 5 

170 008 

2 000 

7. 3 

247 031 

4 000 

9.4 

422 003 

10 000 

12.0 

965 501 

50 000 

16.6 

4 482 024 

100 000 

18. 7 

8 861 127 
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ranging from 20 to 100 000, two minimums 

occur, just as they occur for the previously 

mentioned incompressible case (ref. 5). 

As Q., increases, R first decreases 
1 c 

and then increases monotonically. 

The vertical velocity distribution W 
across the thickness of the layer is illus- 
trated in figure 3. When T is increased, 
the profile is distorted to the left, because 
a thicker layer is less stable, and W 
reaches maximum value sooner and takes a 
longer time to decrease to the vanishing 
value at the upper boundary. 


CONCLUSIONS 


It is found that the compressible 
medium is less stable than the incompres- 
sible medium when heated from below. 

When the rotation rate is small, the inhi- 
bition caused by a magnetic field on the on- 
set of thermal convection is very effective. 
When the rotation rate is large, increasing 
the magnetic field offers only decreasing 
inhibition until a minimum value of the 
critical Rayleigh number is attained. After 
this point, the inhibition increases as the 
magnetic field increases. 



Figure 3. - The z-component velocity W 
profiles for dimensionless thickness 
T = 0. 001 and T = 0. 1. 


Manned Spacecraft Center 

National Aeronautics and Space Administration 
Houston, Texas, November 19, 1970 
188-40-01-01-72 
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APPENDIX 


USE OF APPROXIMATE BOUNDARY CONDITIONS FOR THE 

THIN-LAYER CASE 

For the thin- layer case, the use of approximate boundary conditions can be justi- 
fied from the following comparison in the case without rotation (ref. 3). The exact 
boundary conditions are 

e=0 (Al) 

W = 0 (A2) 

and 

D 2 W + “ DW = 0 at £ = £ 0 , £ 0 + 1 (A3) 

The critical Rayleigh numbers R c are obtained analytically by use of the variational 

method. The critical Rayleigh numbers are then recalculated with the approximate 
boundary conditions 


0=0 

(A4) 

W = 0 

(A5) 

and 


=0 at C = C 0 ,? 0 + 1 

(A6) 

The results are shown in table A-I for Q = 10 and 50. The approximate results differ 
from the exact results by only approximately 2 percent. Hence, the use of approximate 
boundary conditions for the present problem is acceptable. 
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TABLE A-I. - COMPARISON BETWEEN R FROM THE EXACT AND FROM 

c 

THE APPROXIMATE BOUNDARY CONDITIONS 


Q 

■ i 

T 


R 

c 

Error, 

Exact boundary 
conditions 

Approximate 
boundary conditions 

percent 


0. 001 

921. 7 

921. 7 

0 

10 

. 01 

909.4 

909. 3 

-0. 01 


. 1 

803.8 

787. 5 

-2. 00 


0. 001 

1759 

1759 

0 

50 

. 01 

1736 

1736 

0 


. 1 

_ 1 

1535 

1508 

-1. 76 
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